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We use two fundamental theoretical frameworks to study the finite-size (shell) properties of the unitary 
gas in a periodic box: 1) an ab initio Quantum Monte Carlo (qmc) calculation for boxes containing 4 to 
130 particles provides a precise and complete characterization of the finite-size behavior, and 2) a new 
Density Functional Theory (dft) fully encapsulates these effects. The dft predicts vanishing shell 
structure for systems comprising more than 50 particles, and allows us to extrapolate the qmc results 
to the thermodynamic limit, providing the tightest bound to date on the ground-state energy of the 
unitary gas: £,s ^ 0.383(1 ). We also apply the new functional to few-particle harmonically trapped 
systems, comparing with previous calculations. 
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T-he fermion many-body problem plays a fundamen- 
tal role in a vast array of physical systems, from dilute 
gases of cold atoms to nuclear physics in nuclei and 
neutron stars. The universal character of this problem 
- each system is governed by a similar microscopic the- 
ory - coupled with direct experimental access in cold 
atoms, has led to an explosion of recent interest (see 
Refs. [jj for a review). Despite this broad applicability, 
we are far from fully understanding even the simplest 
system: the "unitary gas" comprising equal numbers of 
two fermionic species interacting with a resonant s-wave 
interaction of infinite scattering length q s — > 00. Lack- 
ing any scale beyond the total density n + = n a +n b , 
the unitary gas eschews perturbative expansion and re- 
quires experimental measurement or accurate numerical 
simulation for a quantitative description - the latter is 
presently more precise. Typical Quantum Monte Carlo 
(qmc) calculations, however, can access at most a few 
hundred particles. Density Functional Theory (dft) pro- 
vides a complementary approach through which one 
may extrapolate these results to large systems beyond 
the reach of direct simulation. 

In this Letter, we present the most precise qmc cal- 
culations to date of the unitary gas in a periodic box, 
studying from 4 to 130 particles, thereby providing a 
benchmark for many-body theories. We use this to cali- 
brate a local dft, then use this dft to study the finite-size 
effects ("shell" effects in nuclear physics) and extrapolate 
the qmc results to the thermodynamic limit. We provide 
the most precise bound to date of the universal Bertsch 
parameter [2] £, s = £/£fg ^ 0.383(1). (£ FG = 3/5n + E F 
is the energy density of a free Fermi gas with the same 
total density n+ = kp/(37t 2 ), and Ep = h. 2 k 2 /(2m.) is the 
Fermi energy.) We also explore the finite-size properties 
of the dft - a crucial element in the program to calculate 
properties of finite nuclei with a universal dft [3]. We 
find that a local dft can capture the finite-size effects in 



these systems without the need for particle number pro- 
jection. We limit our discussion to symmetric systems 
(n a = n b ), leaving odd-even staggering to future work, 
as the dft then requires an additional dimensionless 
parameter to characterize the asymmetry n a 7^ n b . 

The qmc results presented here are directly applicable 
to cold 6 Li or 40 K atoms, and constrain dilute neutron 
matter in neutron stars [4]; likewise, the general dft 
approach has myriads of applications throughout cold- 
atom and nuclear physics (see Ref. I5J for a review). 
Our calculation of £,s is consistent with previous results, 
but an order of magnitude more precise. Continuum 
qmc bounds £,$ < 0.40 - 0.44 with an uncertainty no 
better than the last digit lg|l6kio1. Lattice qmc results 
range from « 0.3 - 0.4 IiiHi3jL comparable to ana- 
lytic results [ [14"] . Experimental groups found qualitative 
agreement [15!, which led to precision measurements: 
notably with Duke [16J and Paris [17! quoting 0.39(2) 
and 0.41 (1 ) respectively. 

Dft is an in principle exact approach, widely used 
in in nuclear physics [18J, and in quantum chemistry 
to describe normal (i.e., non-superfluid) systems. It has 
recently been extended to describe the unitary gas |j5j[T9j- 
21]. We build upon one approach - the Superfluid Local 
Density Approximation (slda) - which was originally 
constrained by qmc calculations of the continuum state, 
and then validated with qmc calculations in a harmonic 
trap ri22l[23l (see also Fig.[2j. We focus on translationally 
invariant systems in a periodic box to isolate the finite- 
size effects from the gradient effects. We find that the 
inclusion of an anomalous density is crucial: functionals 
attempting to model the superfluid by adding only gra- 
dient or kinetic corrections [19 , 21] are unable to even 
qualitatively characterize the finite-size effects. 

Our qmc results are based on a fixed-node Diffusion 
Monte Carlo approach that projects out the state of low- 
est energy from the space of all wave functions with 
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Figure i. (color online) Ground-state energy-density £, = 
£/£pg of N+ fermions in a periodic cubic box at the uni- 
tary limit. The circles with error bars are the result of using a 
quadratic least-squares extrapolation to zero effective range of 
our new qmc results. The solid curve is the best fit slda dft. 
The light dotted curve is the functional considered in [19I with 
a = 0.69. For comparison, we have plotted the previous best 
estimate £,s = 0.40(1) (red square) and the current estimate 
£, s = 0.383(1 ) below it to the far right of the figure. Inset: we 
show the typical effective-range dependence £,(k F r e ) with the 
best fit icr error bounds for all-point cubic (solid dark green) 
and five-point quadratic (hatched light yellow) polynomial 
fits. Note that: a) the five-point quadratic model is consistent 
with the full cubic model and has a comparable extrapolation 
error, and b) the inflection point near k F r e « 0.16 necessitates 
a higher-order fit for larger ranges (cubic is sufficient for the 
ranges shown here). Results for N + = 40 show the same 
qualitative behaviour; hence, for the other points we use the 
five-point quadratic extrapolation. 

fixed nodal structure as defined by an initial many -body 
wave function (ansatz). By varying the ansatz, we obtain 
a variational upper bound on the ground-state energy. 
In this work, we use the trial function introduced in (6): 

v' 

where A antisymmetrizes over particles of the same spin 
(either primed or unprimed) and f (r) is a nodeless Jas- 
trow function introduced to reduce the statistical error. 
The antisymmetrized product of s-wave pairing func- 
tions tj)(rtj/) defines the nodal structure: 

4)(r)=^a|| n ||e ik n' r + p(r). (2) 

n 

The sum is truncated (we include ten coefficients) and 
the omitted short-range tail is modelled by the phe- 
nomenological function j3(r) chosen to ensure smooth 
behavior near zero separation. We use the same form 
for p(r) as in |6J with the values b = 0.5 and c = 5, and 
vary the 10 coefficients ocii n n for each IM + to minimize the 
energy as described in Ref. [24]. Representative nodal 
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Table I. Values of £,(N+) = E(N + )/V/£ FG shown in Fig.[i] To 
facilitate comparison with other normalizations in the litera- 
ture, we include the values £,j, 0X (N + ) = E(N + )/E FG (N+) where 
E FG (N + ) is the energy of N+ non-interacting particles in a box. 
The conversion factor is shown in the last column. 

structures are defined by the coefficients in Table [TT] The 
same ansatz suffices for different effective ranges, but 
an independent optimization is required for each N + . 
We simulate the Hamiltonian: 

h 2 ( N+ \ 
K =2 m \-'L V i- 4voVi 2 X sech2 ^ r ij') K (3) 
V k=i i,y / 

with an interspecies interaction of the modified Poschl- 
Teller type (off-resonance intraspecies interactions are 
neglected). We tune to infinite s-wave scattering length 
by setting vo = 1 : the effective range becomes r e = 2/u. 
To extrapolate to the zero-range limit, we simulate at 
u/k F e {12.5,24,36,48,60} for which 0.03 < k F r e < 0.16. 
A careful examination of additional ranges up to k F r e ~ 
0.35 for N + = 40 and N + = 66 (see the inset in Fig. [lj 



N+ q ai a 2 a 3 a 4 a 5 a 6 a 8 a 9 a 10 

10 1600 350 49 16 12 14 14 11 9.0 6.7 

40 160 91 27 0.49 -2.8 -0.086 2.2 2.9 2.5 1.9 

80 -24 13 12 8.2 5.1 3.7 2.7 2.0 1.6 1.0 

120 -51 -17 0.51 7.8 6.3 5.8 4.6 2.5 1.7 1.0 



Table II. Sample coefficients of the pairing function O a||„[| = 
10 _4 ai where I = ||n|| 2 = n 2 +n| + n| = k 2 l 2 /4ni Higher- 
order coefficients are set to zero. 
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reveals that a three-parameter quadratic model in r e 
is necessary and sufficient to extrapolate our results 
without a systematic bias; the results are shown in Fig. [1] 

The energies exhibit definite finite-size effects for 
N + < 50, but are essentially featureless for larger N + . 
This lack of structure is confirmed by the best fit dft (dis- 
cussed below) and disagrees with the results presented 
in Ref. [iol. The values of £, for N + > 50 are distributed 
about the best fit value £,s « 0.383(1), and represent the 
lowest variational bounds to date. Part of the decrease 
from previous results is due to the careful extrapolation 
to zero effective range. The remainder is due to the 
improved optimization of the variational wave function. 

To model the finite-size effects we turn to a local 
dft for the unitary Fermi gas that generalizes the 
Slda originally presented in Ref. li2o1 . In addition to 
the total density n + = 2^ rL |v n | 2 , the slda includes 
both kinetic t + = 2^ n |Vv Tl | 2 and anomalous densities 
v = Xn u Ti v ia- (The + index signifies the sum of the 
contributions coming from the two components a and b; 
u n (r) and v n (r) are the Bogoliubov quasiparticle wave 
functions.) The original three-parameter slda is 

where oc is the inverse effective mass; |3 is the self -energy; 
and y controls the pairing through the regularized cou- 
pling g = 1/(n.y 3 /y — A/a) where A — > oo is a momen- 
tum cutoff that we take to infinity (see Ref. [5] for details). 
One can use numerically the equations for homogeneous 
matter in the thermodynamic to replace the parameters 
p and y with the more physically relevant quantities £,s 
and r| = A/Ep, where A is the pairing gap. 

In principle, the dft can be expressed in terms of 
only the density n + and its gradients. References B21H 
consider local formulations of this type (called Ex- 
tended Thomas-Fermi (etf) functionals). Since gradients 
vanish in the periodic box, etf functionals reduce to 
£etf(ti+) = £,s£fg an d exhibit no finite-size structure, 
contrary to the Qmc results. Reference [19] adds oct+, 
but without v^v the finite-size effects do not correlate 
with the qmc behavior (see Fig. [lj and the best fit to 
our results is also flat (oc — > 0). Furthermore, such func- 
tionals cannot qualitatively reproduce the quasiparticle 
dispersion relationship, an attractive feature of the slda 
(see also Ref. ps\). 

The best fit three-parameter slda functional Q - a = 
1.26(2), £, s = 0.3826(5), and n = 0.87(2) - is shown in 
Fig. [1] It fits the 23 qmc points from N + = 4toN + = 
130 with a reduced chi squared Xr = 0.7, indicating 
complete consistency. Although remarkable, the fit is 
not completely satisfactory: 1) Fitting the exact two- 
particle energy £,2 = — 0.41 53 - ■ raises Xr = 2-0, and 
2) the best fit gap parameter n and inverse effective 
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Figure 2. (color online) Ground-state energy of the harmon- 
ically trapped unitary Fermi gas (in units where ftu> = 1) 
scaled to demonstrate the asymptotic form 16E 2 /(3N + ) 8 / 3 = 
£,s(1 + cx + 0(x 7 / 6 )) 2 predicted by the low-energy effective 
theory of Ref. [32) . The best fit slda (solid blue line) is com- 
pared with zero-range results for N + e {4,6} from Ref. [33] , 
and finite-range qmc results from Ref. [23! (upper red dots) 
and Ref. [34I (green pluses). The latter have significantly lower 
energy despite having a slightly large effective range, sug- 
gesting that the wave functions in Ref. [23I were not fully 
optimized. We expect careful optimization and zero-range 
extrapolation to bring the qmc results for large N + in line with 
the dft as discussed in the text. 



mass a are inconsistent with the values r\ = 0.50(5) and 
a = 1 .09(2) obtained from the N + = 66 qmc quasiparticle 
dispersion relation [8 26J, and the values n = 0.45(5) [27] 
and n = 0.44(3) [28J extracted from experimental data. 

These deficiencies might be remedied by generalizing 
the slda. As noted in Ref. [5], the following combination 
of divergent kinetic and anomalous densities is finite: 



h 2 T + q + 



h 2 T+ 



y^y 



otn/ 3 /y-A 



2m a 2m 
The lack of scales thus dictates the functional form: 



(5) 



£(K,u+) = UQ) £ FG (n+), Q = K/£ FG (n+), (6) 

where Q is dimensionless, and the regularization con- 
dition depends on Q through the function y(Q)- The 
original slda is linear f,(Q) = ocQ + (3 with constant 
y(Q) — Y- This generalized functional can fit any mono- 
tonic £,(N + ), including the exact N + = 2 point. For 
N + > 6, £,(N + ) is not monotonic and the functional is 
in principle constrained. For example, requiring that 
£, = £,s at both N + w 6.2(2) and N + = 00 fixes the ratio 
r\/oc = 0.69(2). (As an aside, we note that the momentum 
distribution in the dft relates this to the "contact" C: 
ri/ct = VlC/^l ~ 0-44 - 0.49; see Refs. [29] and references 
therein, though it is not clear that this property should 
be trusted.) In practice, the errors and the discreteness 
in N + leave room for flexibility in the functional form, 
and we have found several generalized functional forms 
with Xr ~ 1 -5 while constraining r\ = 0.50. We may have 
to accept the discrepancy in a as a limitation of the dft. 

However, generalizing the slda may not be needed: 
analyzing the "symmetric heavy-light ansatz" [30], (jus- 
tified by lattice qmc calculations B12I ), we find that the 
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simple three-parameter slda suffices (x? ~ 0-5) with 
reasonable a = 0.96(2), r, = 0.51(1), and £, = 0.322(2) - 
slightly higher than the £, = 0.31 (1 ) extracted in [30!. 

It is not trivial that the simple dft Q captures all 
finite-size effects above N + = 4 to high precision in both 
calculations, indicating that the slda may be used to 
extrapolate to the thermodynamic limit. We note that 
no particle-number projection is required - a quite ill- 
defined procedure often considered necessary in nuclear 
physics [ j3~rj : Perhaps improved nuclear functionals may 
similarly capture finite-size effects through local anoma- 
lous densities in the spirit of v. 

To finish, we consider harmonically trapped systems 
in Fig. [2] As discussed in [32], the energy may be ex- 
pressed asE(N + ) = lhci) V / ^(3N + ) 4 / 3 (1 +cx + 0(x 7 / 6 )) 
where x = (3N + )~ 2 / 3 and c is expressed in terms of low- 
energy coefficients. As demonstrated by the zero-range 
N + 6 {4,6} results of [33], the dft still over-estimates the 
energy for small systems, most likely because we have 
omitted the gradient terms in the functional that vanish 
in homogeneous systems. 

For large N + the dft has the expected asymptotic form 
with intercept £,s = 0.383 unlike the finite-range QMC 
results of Refs. f23~}|34"j . This is qualitatively consistent 
with the leading effective-range corrections which scale 
asymptotically as x -1 / 4 ; the systematic overestimation 
of the energy by the variational qmc approach might also 
contribute. We defer further discussion until carefully 
extrapolated zero-range results are published. 

To summarize, we present the most precise Quantum 
Monte Carlo calculations to date of a symmetric unitary 
Fermi gas in a periodic box comprising 4 to 130 particles. 
By carefully characterizing and extrapolating these re- 
sults to zero effective range, we have completely mapped 
out the finite-size effects. These results are used to ana- 
lyze the structure of a Density Functional Theory for the 
symmetric unitary gas, and it is shown that the simplest 
three-parameter form of Eq. Q fully accounts for all 
shell effects to within the statistical errors of the qmc 
results without the need for particle-number projection; 
a more complicated form, however, may be required 
to capture both the finite-size effects and the quasipar- 
ticle dispersions. The dft predicts no significant shell 
corrections beyond 50 particles, and the qmc calcula- 
tions confirm this, allowing us to extract a precise upper 
bound on the universal equation of state £,s ^ 0.383(1 ), 
an order of magnitude improvement in precision over 
previous bounds and the lowest bound of any varia- 
tional method to date. The functional in its latest form 
is well constrained, but leads to slight disagreements 
with qmc predictions for harmonic traps. Converging 
both qmc and dft approaches promises to be a fruitful 
direction of future research. 
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